#ifndef _BER_H_
#define _BER_H_

#include <vector>
#include <complex>
#include <algorithm>
#include <functional>
#include <numeric>
#include <iostream>
#include "complex_ops.h"
#include "criterion.h"
#include "NumericTypes.h"

template
<
    typename SignalToNoiseRate = Double<1, 0>,
    long AttemptCount = 1000
>
class BitErrorRate
{
public:

    const char* getName()
    {
        const static char name[] = "BitErrorRate"; 
        return name;
    }
    
    template
    <
        typename SignalSystem
    >
    double operator()(SignalSystem& system,
        double snr = SignalToNoiseRate().value,
        double count = AttemptCount)
    {
        using namespace std;
        
        if(system.empty())
            return 0;
        
        vector<complex<double> > noise(system[0].size());
    
        double errors = 0;
    
        for (int i = 0; i < system.size(); ++ i)
        { 
            for (int j = 0; j < count / system.size(); ++ j)
            {
                auto recv = system[i];
            
                transform(
                    system[i].begin(), system[i].end(),
                    recv.begin(),
                    bind2nd(plus_random(), snr)
                );
            
                // system
                double dista = BitErrorRate::dist(system[0], recv);
                int opt = 0;
                for(int p = 1; p < system.size(); ++ p)
                {
                    double newDist = BitErrorRate::dist(system[p], recv);
                    if(newDist < dista)
                    {
                        opt = p;
                        dista = newDist;
                    }
                }
            
                if(opt != i)
                    ++ errors;
            }
        }
    
        return errors / count; 
    }
    
private:
    
    static double dist(
        std::vector<std::complex<double> > x,
        std::vector<std::complex<double> > y
    )
    {
        using namespace std;
        
        if(x.size() > y.size())
            swap(x, y);
        
        auto diff = y;
        
        transform(
            x.begin(), x.end(),
            diff.begin(), diff.begin(),
            minus<complex<double> >()
        );
        
        vector<double> rslt(diff.size());

        transform(
            diff.begin(), diff.end(),
            rslt.begin(),
            absolute<complex<double>,double>()
        );

        std::binder2nd
        <
            power<int, double>
        >
        power2(
            power<int, double>(),
            2
        );
        
        transform(
            rslt.begin(), rslt.end(),
            rslt.begin(),
            power2
        );
      
        double sum =
            accumulate(
                rslt.begin(), rslt.end(),
                0.0
            );
        
        return sqrt(sum);
    }
};

#endif // _BER_H_
